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Glossary 

• Time series: One dimensional array of numbers {xi), i = 1,...,N, 



representing values of an observable x usually measured equidistant (or 
nearly equidistant) in time. 

• Complex system: A system consisting of many non-linearly interact- 
ing components. It cannot be split into simpler sub-systems without 
tampering the dynamical properties. 

• Scaling law: A power law with a scaling exponent (e. g. a) describing 
the behaviour of a quantity F (e. g., fluctuation, spectral power) as 
function of a scale parameter s (e. g., time scale, frequency) at least 
asymptotically: F{s) ~ s". The power law should be valid for a large 
range of s values, e. g., at least for one order of magnitude. 

• Fractal system: A system characterised by a scaling law with a frac- 
tal, i. e., non-integer exponent. Fractal systems are characterised by 
self-similarity, i. e., a magnification of a small part is statstically equiv- 
alent to the whole. 
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• Self-aifRne system: Generalization of a fractal system, where different 
magnifications s and s' = have to be used for different directions in 
order to obtain a statistically equivalent magnihcation. The exponent 
H is called Hurst exponent. Self-affine time series and time series 
becoming self-affine upon integration are commonly denoted as fractal 
using a less strict terminology 

• Multifractal system: A system characterised by scaling laws with an 
infinite number of different fractal exponents. The scaling laws must 
be valid for the same range of the scale parameter. 

• Crossover: Change point in a scaling law, where one scaling exponent 
applies for small scale parameters and another scaling exponent applies 
for large scale parameters. The center of the crossover is denoted by 
its characteristic scale parameter Sx in this article. 

• Persistence: In a persistent time series, a large value is usually (i. e., 
with high statistical preference) followed by a large value and a small 
value is followed by a small value. A fractal scaling law holds at least 
for a limited range of scales. 

• Short-term correlations: Correlations that decay sufficiently fast 
that they can be described by a characteristic correlation time scale; 
e. g., exponentially decaying correlations. A crossover to uncorrelated 
behaviour is observed on larger scales. 

• Long-term correlations: Correlations that decay sufficiently slow 
that a characteristic correlation time scale cannot be defined; e. g., 
power-law correlations with an exponent between and 1. Power-law 
scaling is observed on large time scales and asymptotically. The term 
long-range correlations should be used if the data is not a time series. 

• Non-stationarities: If the mean or the standard deviation of the data 
values change with time, the weak definition of stationarity is violated. 
The strong definition of stationarity requires that all moments remain 
constant, i. e., the distribution density of the values does not change 
with time. Non-stationarities like monotonous, periodic, or step-like 
trends are often caused by external effects. In a more general sense, 
changes in the dynamics of the system also represent non-stationarities. 
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1 Definition of the Subject and Its Impor- 
tance 



Data series generated by complex systems exhibit fluctuations on a wide 
range of time scales and/or broad distributions of the values. In both equi- 
librium and non-equilibrium situations, the natural fluctuations are often 
found to follow a scaling relation over several orders of magnitude. Such 
scaling laws allow for a characterisation of the data and the generating com- 
plex system by fractal (or multifractal) scaling exponents, which can serve as 
characteristic fingerprints of the systems in comparisons with other systems 
and with models. Fractal scaling behaviour has been observed, e. g., in many 
data series from experimental physics, geophysics, medicine, physiology, and 
even social sciences. Although the underlying causes of the observed fractal 
scaling are often not known in detail, the fractal or multifractal characteri- 
sation can be used for generating surrogate (test) data, modelling the time 
series, and deriving predictions regarding extreme events or future behaviour. 
The main application, however, is still the characterisation of different states 
or phases of the complex system based on the observed scaling behaviour. 
For example, the health status and different physiological states of the hu- 
man cardiovascular system are represented by the fractal scaling behaviour 
of the time series of intervals between successive heartbeats, and the coars- 
ening dynamics in metal alloys are represented by the fractal scaling of the 
time-dependent speckle intensities observed in coherent X-ray spectroscopy. 

In order to observe fractal and multifractal scaling behaviour in time 
series, several tools have been developed. Besides older techniques assum- 
ing stationary data, there are more recently established methods differen- 
tiating truly fractal dynamics from fake scaling behaviour caused by non- 
stationarities in the data. In addition, short-term and long-term correlations 
have to be clearly distinguished to show fractal scaling behaviour Tinam- 
biguously. This article describes several methods originating from Statistical 
Physics and Applied Mathematics, which have been used for fractal and 
multifractal time series analysis in stationary and non-stationary data. 

2 Introduction 

The characterisation and understanding of complex systems is a difficult task, 
since they cannot be split into simpler subsystems without tampering the dy- 
namical properties. One approach in studying such systems is the recording 
of long time series of several selected variables (observables) , which refiect 
the state of the system in a dimensionally reduced representation. Some 
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systems are characterised by periodic or nearly periodic behaviour, which 
might be caused by oscillatory components or closed-loop regulation chains. 
However, in truly complex systems such periodic components are usually not 
limited to one or two characteristic frequencies or frequency bands. They 
rather extend over a wide spectrum, and fluctuations on many time scales as 
well as broad distributions of the values are found. Often no specific lower 
frequency limit - or, equivalently, upper characteristic time scale - can be 
observed. In these cases, the dynamics can be characterised by scaling laws 
which are valid over a wide (possibly even unlimited) range of time scales 
or frequencies; at least over orders of magnitude. Such dynamics are usu- 
ally denoted as fractal or multifractal, depending on the question if they are 
characterised by one scaling exponent or by a multitude of scaling exponents. 

The first scientist who applied fractal analysis to natural time series is 
Benoit B. Mandelbrot [HEIE], who included early approaches by H.E. Hurst 
regarding hydrological systems [HE]. For extensive introductions describing 
fractal scaling in complex systems, we refer to |6l El [U |9l [TOl HH [111 |13]. In 
the last decade, fractal and multifractal scaling behaviour has been reported 
in many natural time series generated by complex systems, including 

• geophysics time series (recordings of temperature, precipitation, water 
runoff, ozone levels, wind speed, seismic events, vegetational patterns, 
and climate dynamics), 

• medical and physiological time series (recordings of heartbeat, respi- 
ration, blood pressure, blood flow, nerve spike intervals, human gait, 
glucose levels, and gene expression data), 

• DNA sequences (they are not actually time series) , 

• astrophysical time series (X-ray light sources and sunspot numbers), 

• technical time series (internet traffic, highway traffic, and neutronic 
power from a reactor), 

• social time series (finance and economy, language characteristics, fatal- 
ities in conflicts), as well as 

• physics data (also going beyond time series), e. g., surface roughness, 
chaotic spectra of atoms, and photon correlation spectroscopy record- 
ings. 

If one flnds that a complex system is characterised by fractal (or multi- 
fractal) dynamics with particular scaling exponents, this flnding will help in 
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obtaining predictions on the future behaviour of the system and on its reac- 
tion to external perturbations or changes in the boundary conditions. Phase 
transitions in the regulation behaviour of a complex system are often associ- 
ated with changes in their fractal dynamics, allowing for a detection of such 
transitions (or the corresponding states) by fractal analysis. One example for 
a successful application of this approach is the human cardiovascular system, 
where the fractality of heartbeat interval time series was shown to reflect 
certain cardiac impairments as well as sleep stages [HI |T3] . In addition, one 
can test and iteratively improve models of the system until they reproduce 
the observed scaling behaviour. One example for such an approach is climate 
modelling, where the models were shown to need input from volcanos and so- 
lar radiation in order to reproduce the long-term correlated (fractal) scaling 
behaviour (TU] previously found in observational temperature data [T7] . 

Fractal (or multifractal) scaling behaviour certainly cannot be assumed a 
priori, but has to be established. Hence, there is a need for refined analysis 
techniques, which help to differentiate truly fractal dynamics from fake scal- 
ing behaviour caused, e. g., by non-stationarities in the data. If conventional 
statistical methods are applied for the analysis of time series representing 
the dynamics of a complex system \W\ . there are two major problems, 
(i) The number of data series and their durations (lengths) are usually very 
limited, making it difficult to extract significant information on the dynam- 
ics of the system in a reliable way. (ii) If the length of the data is extended 
using computer-based recording techniques or historical (proxy) data, non- 
stationarities in the signals tend to be superimposed upon the intrinsic fluc- 
tuation properties and measurement noise. Non-stationarities are caused by 
external or internal effects that lead to either continuous or sudden changes 
in the average values, standard deviations or regulation mechanism. They 
are a major problem for the characterisation of the dynamics, in particular 
for finding the scaling properties of given data. 

Following the description of important properties of fractal and multi- 
fractal time series and the definition of our primary quantities of interest in 
Section 3, we focus on methods for the analysis of self-affine or (mono-)fractal 
data in Sections 4 and 5. While Section 4 describes traditional approaches 
for the analysis of stationary time series. Section 5 discusses more recently 
established methods applicable for non- stationary data. Section 6 describes 
techniques for multifractal time series analysis. The consequences of fractal 
scaling behaviour on the statistics of extreme events and their return inter- 
vals are presented in Section 7, and a few standard models for fractal and 
multifractal time series are described in Section 8 before an outlook on future 
directions in Section 9. 
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3 Fractal and Multi fractal Time Series 



3.1 Fractality, Self- Affinity, and Scaling 

The topic of this article is the fractahty (and/or multifractahty) of time 
series. Since fractals and multifractals in general are discussed in many other 
articles of the encyclopedia, the concept is not thoroughly explained here. In 
particular, we refer to the articles ... and ... for the formalism describing 
fractal and multifractal structures, respectively. 

In a strict sense, most time series are one dimensional, since the values 
of the considered observable are measured in homogeneous time intervals. 
Hence, unless there are missing values, the fractal dimension of the support 
is -D(O) = 1. However, there are rare cases where most of the values of a 
time series are very small or even zero, causing a dimension D(0) < 1 of 
the support. In these cases, one has to be very careful in selecting appropri- 
ate analysis techniques, since many of the methods presented in this article 
are not accurate for such data; the Wavelet Transform Modulus Maxima 
technique (see Section 6.2) is the most advanced applicable method. 

Even if the fractal dimension of support is one, the information dimension 
D{1) and the correlation dimension D(2) can be studied. As we will see 
in Section 6.1, D{2) is in fact explicitly related to all exponents studied 
in monofractal time series analysis. However, usually a slightly different 
approach is employed based on the notion of self-affinity instead of (multi-) 
fractality. Here, one takes into account that the time axis and the axis of 
the measured values x{t) are not equivalent. Hence, a rescaling of time t by 
a factor a may require rescaling of the series values x{t) by a different factor 

in order to obtain a statistically similar (i. e., self-similar) picture. In this 
case the scaling relation 



holds for an arbitrary factor a, describing the data as self-affine (see, e. g., [B]). 
The i/ttrsi exponent H (after the water engineer H.E. Hurst |1]) characterises 
the type of self affinity. Figure 1(a) shows several examples of self-affine time 
series with different H. The trace of a random walk (Brownian motion, third 
line in Fig. 1(a)), for example, is characterised by if = 0.5, implying that 
the position axis must be rescaled by a factor of 2 if the time axis is rescaled 
by a factor of 4. Note that self-affine series are often denoted as fractal 
even though they are not fractal in the strict sense. In this article the term 
"fractal" will be used in the more general sense including all data, where a 
Hurst exponent H can be reasonably defined. 

The scaling behaviour of self-affine data can also be characterised by look- 
ing at their mean-square displacement. Since the mean-square displacement 
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of a random walker is known to increase linear in time, {x^{t)) ~ t, deviations 
from this law will indicate the presence of self-affine scaling. As we will see 
in Section 4.4, one can thus retrieve the Hurst (or self-affinity) exponent H 
by studying the scaling behaviour of the mean-square dispalcement, or the 
mean-square fluctuations {x^{t)) ~ t^^ . 

3.2 Persistence, Long- and Short-term Correlations 

Self-affine data are persistent in the sense that a large value is usually (i. e., 
with high statistical preference) followed by a large value and a small value 
is followed by a small value. For the trace of a random walk, persistence 
on all time scales is trivial, since a later position is just a former one plus 
some random increment (s). The persistence holds for all time scales, where 
the self-affinity relation ([!]) holds. However, the degree of persistence can 
also vary on different time scales. Weather is a typical example: while the 
weather tomorrow or in one week is probably similar to the weather today 
(due to a stable general weather condition), persistence is much harder to be 
seen on longer time scales. 

Considering the increments Axj = Xj — of a self-affine series, (xj), 
i = 1, . . . ,N with N values measured equidistant in time, one finds that the 
Axj can be either persistent, independent, or anti-persistent. Examples for 
all cases are shown in Fig. 1(b). In our example of the random walk with 
H = 0.5 (third line in the figure), the increments (steps) are fully independent 
of each other. Persistent and anti-persistent increments, where a positive 
increment is likely to be followed by another positive or negative increment, 
respectively, are also leading to persistent integrated series Xi = Z]}=i ^^j- 

For stationary data with constant mean and standard deviation the auto- 
covariance function of the increments, 

^ N-s 

C{s) = {Axi Axi+,) = — — ^""i+s- (2) 

^ i=i 

can be studied to determine the degree of persistence. If C(s) is divided 
by the variance ((Axj)^), it becomes the auto-correlation function; both are 
identical if the data are normalised with unit variance. If the Axi are un- 
correlated (as for the random walk), C{s) is zero for s > 0. Short-range 
correlations of the increments Axj are usually described by C{s) declining 
exponentially, 

C(s) ~exp(-s/tx) (3) 
with a characteristic decay time tx- Such behaviour is typical for increments 
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generated by an auto-regressive (AR) process 



Axj = cAa;j_i + £j 



(4) 



with random uncorrelated offsets £j and c = exp(— l/^x)- Figure 2(a) sfiows 
tlie auto-correlation function for one configuration of an AR process with 



For so-called long-range correlations J^C{s)ds diverges in the limit of 
infinitely long series {N — > oo). In practice, this means that cannot be 
defined because it increases with increasing N. For example, C(s) declines 
as a power-law 



with an exponent < 7 < 1. Figure 2(b) shows C{s) for one configuration 
with 7 = 0.4. This type of behaviour can be modelled by the Fourier filtering 
technique (see Section 8.1). Long-term correlated, i. e. persistent, behaviour 
of the Axi leads to self-affine scahng behaviour of the Xi, characterised by 
= 7/2, as will be shown below. 

3.3 Crossovers and Non-stationarities in Time Series 

Short-term correlated increments Axi characterised by a finite characteristic 
correlation decay time tx lead to a crossover in the scahng behaviour of the 
integrated series ,Xj = Yl]=iAxj, see Fig. 2(a) for an example. Since the 
position of the crossover might be numerically different from t-^, we denote 
it by Sx here. Time series with a crossover are not self-affine and there is no 
unique Hurst exponent H characterising them. While H > 0.5 is observed on 
small time scales (indicating correlations in the increments) , the asymptotic 
behaviour (for large time scales s ':$> and 3> Sx) is always characterised 
hj H = 0.5, since all correlations have decayed. Many natural recordings are 
characterised by pronounced short-term correlations in addition to scaling 
long-term correlations. For example, there are short-term correlations due 
to particular general weather situations in temperature data and due to res- 
pirational effects in heartbeat data. Crossovers in the scaling behaviour of 
complex time series can also be caused by different regulation mechanisms 
on fast and slow time scales. Fluctuations of river runoff, for example, show 
different scaling behaviour on time scales below and above approximately 
one year. 

Non-stationarities can also cause crossovers in the scaling behaviour of 
data if they are not properly taken into account. In the most strict sense, 
non-stationarities are variations in the mean or the standard deviation of 
the data (violating weak stationarity) or the distribution of the data values 



ix = 48. 
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(violating strong stationarity). Non-stationarities like monotonous, periodic 
or step- like trends are often caused by external effects, e. g., by the greenhouse 
warming and seasonal variations for temperature records, different levels of 
activity in long-term physiological data, or unstable light sources in photon 
correlation spectroscopy. Another example for non-stationary data is a record 
consisting of segments with strong fluctuations alternating with segments 
with weak fluctuations. Such behaviour will cause a crossover in scaling 
at the time scale corresponding to the typical duration of the homogeneous 
segments. Different mechanisms of regulation during different time segments 
- like, e. g., different heartbeat regulation during different sleep stages at 
night - can also cause crossovers; they are regarded as non-stationarities 
here, too. Hence, if crossovers in the scaling behaviour of data are observed, 
more detailed studies are needed to find out the cause of the crossovers. 
One can try to obtain homogenous data by splitting the original series and 
employing methods that are at least insensitive to monotonous (polynomially 
shaped) trends. 

To characterise a complex system based on time series, trends and fluctu- 
ations are usually studied separately (see, e. g., |2D] for a discussion). Strong 
trends in data can lead to a false detection of long-range statistical persis- 
tence if only one (non-detrending) method is used or if the results are not 
carefully interpreted. Using several advanced techniques of scaling time se- 
ries analysis (as described in Chapter 5) crossovers due to trends can be 
distinguished from crossovers due to different regulation mechanisms on fast 
and slow time scales. The techniques can thus assists in gaining insight into 
the scaling behaviour of the natural variability as well as into the kind of 
trends of the considered time series. 

It has to be stressed that crossovers in scaling behaviour must not be 
confused with multifractality. Even though several scaling exponents are 
needed, they are not applicable for the same regime (i. e., the same range 
of time scales). Real multifractality, on the other hand, is characterised by 
different scaling behaviour of different moments over the full range of time 
scales (see next section). 

3.4 Multifractal Time Series 

Many records do not exhibit a simple monofractal scaling behaviour, which 
can be accounted for by a single scaling exponent. As discussed in the previ- 
ous section, there might exist crossover (time-) scales Sx separating regimes 
with different scaling exponents. In other cases, the scaling behaviour is 
more complicated, and different scaling exponents are required for different 
parts of the series. In even more complicated cases, such different scaling 
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behaviour can be observed for many interwoven fractal subsets of the time 
series. In this multitude of scaling exponents is required for a full 

description of the scaling behaviour in the same range of time scales, and a 
multifractal analysis must be applied. 

Two general types of multifractality in time series can be distinguished: 
(i) Multifractality due to a broad probabihty distribution (density function) 
for the values of the time series, e. g. a Levy distribution. In this case the 
multifractality cannot be removed by shuffling the series, (ii) Multifractality 
due to different long-term correlations of the small and large fluctuations. 
In this case the probability density function of the values can be a regu- 
lar distribution with flnite moments, e. g., a Gaussian distribution. The 
corresponding shuffled series will exhibit non- multifractal scaling, since all 
long-range correlations are destroyed by the shuffling procedure. Randomly 
shuffling the order of the values in the time series is the easiest way of gen- 
erating surrogate data; however, there are more advanced alternatives (see 
Chapter 8). If both kinds of multifractality are present, the shuffled series 
will show weaker multifractality than the original series. 

A multifractal analysis of time series will also reveal higher order cor- 
relations. Multifractal scaling can be observed if, e. g., three or four-point 
correlations scale differently from the standard two-point correlations stud- 
ied by classical autocorrelation analysis (Eq. (|2|). In addition, multifractal 
scaling is observed if the scaling behaviour of small and large fluctuations is 
different. For example, extreme events might be more or less correlated than 
typical events. 

4 Methods for Stationary Fractal Time Se- 
ries Analysis 

In this chapter we describe four traditional approaches for the fractal anal- 
ysis of stationary time series, see [211 E2l EH] for comparative studies. The 
main focus is on the determination of the scaling exponents H or 7, defined in 
Eqs. ([1]) and (|5]), respectively, and linked byif = l— 7/2in long-term persis- 
tent data. Methods taking non-stationarities into account will be discussed 
in the next chapter. 

4.1 Autocorrelation Function Analysis 

We consider a record (xj) oii = 1, . . . , N equidistant measurements. In most 
applications, the index i will correspond to the time of the measurements. 
We are interested in the correlation of the values Xi and Xi+s for different 
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time lags, i. e. correlations over different time scales s. In order to remove a 
constant offset in the data, the mean (x) = J2iLi is usually subtracted, 
Xi = Xi — (x). Alternatively, the correlation properties of increments Xi = 
Axi = Xi — of the original series can be studied (see also Section 3.2). 
Quantitatively, correlations between x-values separated by s steps are defined 
by the (auto-) covariance function C{s) = {xiXi+s) or the (auto-) correlation 
function C{s)/{xf), see also Eq. 

As already mentioned in Section 3.2, the short-term correlated if 

C{s) declines exponentially, C{s) ~ exp(— s/tx), and long-term correlated if 
C{s) declines as a power-law C{s) oc s~~^ with a correlation exponent < 
7 < 1 (see Eqs. ^ and ([s]), respectively). As illustrated by the two examples 
shown in Fig. 2, a direct calculation of C{s) is usually not appropriate due to 
noise superimposed on the data Xi and due to underlying non-stationarities of 
unknown origin. Non-stationarities make the definition of C{s) problematic, 
because the average (x) is not well-defined. Furthermore, C{s) strongly 
fluctuates around zero on large scales s (see Fig. 2(b)), making it impossible 
to find the correct correlation exponent 7. Thus, one has to determine the 
value of 7 indirectly. 

4.2 Spectral Analysis 

If the time series is stationary, we can apply standard spectral analysis tech- 
niques (Fourier transform) and calculate the power spectrum S{f) of the time 
series {xi) as a function of the frequency / to determine self-affine scaling be- 
haviour [23]. For long-term correlated data characterised by the correlation 
exponent 7, we have 

Sif)-r^ with /3 = l-7. (6) 

The spectral exponent f3 and the correlation exponent 7 can thus be obtained 
by fitting a power-law to a double logarithmic plot of the power spectrum 
S{f). An example is shown in Fig. 3. The relation (|6| can be derived from 
the Wiener-Khinchin theorem (see, e. g., [25]). If, instead of Xi = Axi the 
integrated runoff time series is Fourier transformed, i. e., Xi = XiY^'^j^i Axj, 
the resulting power spectrum scales as S{f) ~ f~'^~^ ■ 

Spectral analysis, however, does not yield more reliable results than auto- 
correlation analysis unless a logarithmic binning procedure is applied to the 
double logarithmic plot of S{f) [21j, see also Fig. 3. I. e., the average of 
log5'(/) is calculated in successive, logarithmically wide bands from a"/o 
to a"+^/o, where /o is the minimum frequency, a > 1 is a factor (e. g., 
a = 1.1), and the index n is counting the bins. Spectral analysis also requires 
stationarity of the data. 



12 



4.3 Hurst's Rescaled-Range Analysis 

The first method for the analysis of long-term persistence in time series based 
on random walk theory has been proposed by the water construction engineer 
Harold Edwin Hurst (1880-1978), who developed it while working in Egypt. 
His so-called rescaled range analysis {R/S analysis) [H E], IH El [6] begins 
with splitting of the time series (x,) into non-overlapping segments u of size 
(time scale) s (first step), yielding A^^ = int(A^/s) segments altogether. In 
the second step, the profile (integrated data) is calculated in each segment 
= 0, . . . , AT, - 1, 

j 3 ' s 

"^(j) ^ ^ (-^i/s+j ip^us+i) s) ^ ^ •^vs+i ~ ^ ^ -^us+i' (7) 
i=l i=l i=l 

By the subtraction of the local averages, piecewise constant trends in the 
data are eliminated. In the third step, the differences between minimum and 
maximum value {ranges) Ry{s) and the standard deviations Sj,{s) in each 
segment are calculated. 



\ 



Finally, the rescaled range is averaged over all segments to obtain the fluc- 
tuation function F{s), 

Fnsis) = ^ ^ - for . » 1, (9) 

where H is the Hurst exponent already introduced in Eq. ([T]). One can show 
[B [21] that H is related to /3 and 7 by 2i/ ^ 1 + /? = 2 - 7 (see also Eqs. (g 



and (14)). Note that < 7 < 1, so that the right part of the equation does 
not hold unless 0.5 < H < 1. The relationship does not hold in general for 
multifractal data. Note also that H actually characterises the self-affinity of 
the profile function ([T]), while /3 and 7 refer to the original data. 

The values of H, that can be obtained by Hurst's rescaled range analysis, 
are limited to < if < 2, and significant inaccuracies are to be expected 
close to the bounds. Since H can be increased or decreased by 1 if the data is 
integrated {xj J2i=iS:i) or differentiated (xj Xi — Xi-i), respectively, one 
can always find a way to calculate H by rescaled range analysis provided the 
data is stationary. While values H < 1/2 indicate long-term anti-correlated 
behaviour of the data Xj, H > 1/2 indicates long-term positively correlated 
behaviour. For power-law correlations decaying faster than 1/s, we have 
if = 1/2 for large s values, like for uncorrelated data. 
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Compared with spectral analysis, Hurst's rescaled range analysis yields 
smoother curves with less effort (no binning procedure is necessary) and 
works also for data with piecewise constant trends. 



4.4 Fluctuation Analysis (FA) 

The standard fluctuation analysis (FA) [HI |26] is also based on random walk 
theory. For a time series (xi), i = 1, . . . , N, with zero mean, we consider the 
global profile, i. e., the cumulative sum (cf. Eq. ([T])) 



yiJ) = ^2^^^ j = o,i,2,...,iv, (10) 



i=l 



and study how the fluctuations of the proflle, in a given time window of size 
s, increase with s. The procedure is illustrated in Fig. 4 for two values of 
s. We can consider the proflle Y{j) as the position of a random walker on 
a linear chain after j steps. The random walker starts at the origin and 
performs, in the zth step, a jump of length Xi to the right, if Xj is positive, 
and to the left, if Xi is negative. 

To flnd how the square-fluctuations of the proflle scale with s, we flrst 
divide each record of elements into Ng = int(A^/s) non-overlapping seg- 
ments of size s starting from the beginning (see Fig. 4) and another Ng 
non-overlapping segments of size s starting from the end of the considered 
series. This way neither data at the end nor at the beginning of the record 
is neglected. Then we determine the fluctuations in each segment u. 

In the standard FA, we obtain the fluctuations just from the values of the 
proflle at both endpoints of each segment u = 0, . . . , Ng — 1, 

F2^(z.,.) = [F(z..)-F((z/ + l)s)]^ (11) 

(see Fig. 4) and analogous for u = Ng, . . . , 2Ng — 1, 

F^^{u, s) = [Y{N -{u- Ng)s) -Y{N-{iy + l- Ng)s)]\ (12) 

Then we average Fp^^lu, s) over all subsequences to obtain the mean fluctu- 
ation F2{s), 



Fois) 



1 2Ns-l 



,1/2 



s". (13) 



By deflnition, F2{s) can be viewed as the root-mean-square displacement of 
the random walker on the chain, after s steps (the reason for the index 2 will 
become clear later). For uncorrected Xi values, we obtain Fick's diffusion 
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law F2{s) ~ s^/^. For the relevant case of long-term correlations, where C{s) 
follows the power-law behaviour of Eq. (|5]), F2{s) increases by a power law, 

F2{s) ~ with a^H, (14) 

where the fluctuation exponent a is identical with the Hurst exponent H for 
mono-fractal data and related to 7 and /3 by 

2a = l + /5 = 2-7. (15) 

The typical behaviour of -^2(5) for short-term correlated and long-term cor- 



related data is illustrated in Fig. 2. The relation (15) can be derived straight- 
forwardly by inserting Eqs. (10), ([2]), and (|5| into Eq. (11) and separating 
sums over products XiXj with identical and different i and j, respectively. 

The range of the a values that can be studied by standard FA is limited to 
< a < 1, again with significant inaccuracies close to the bounds. Regarding 
integration or differentiation of the data, the same rules apply as listed for 
H in the previous section. The results of FA become statistically unreliable 
for scales s larger than one tenth of the length of the data, i. e. the analysis 
should be limited by s < A^/10. 



5 Methods for Non- Stationary Fractal Time- 
Series Analysis 

5.1 Wavelet Analysis 

The origins of wavelet analysis come from signal theory, where frequency de- 
compositions of time series were studied [271 [2S]- Like the Fourier transform, 
the wavelet transform of a signal x{t) is a convolution integral to be replace 
by a summation in case of a discrete time series {xi),i = 1, . . . , N , 

I foo I ^ 

L^{t,s) = - x{t)tP[{t-r)/s]dt = -J2xi^[{t-T)/s]. (16) 

S J— 00 S 

Here, ip{t) is a so-called mother wavelet, from which all daughter wavelets 
ipT,s{t) = 'ip{{t — t)/s evolve by shifting and stretching of the time axis. The 
wavelet coefficients L^{t, s) thus depend on both, time position r and scale 
s. Hence, the local frequency decomposition of the signal is described with 
a time resolution appropriate for the considered frequency / = 1/s (i. e., 
inverse time scale). 
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All wavelets tp{t) must have zero mean. They are often chosen to be 
orthogonal to polynomial trends, so that the analysis method becomes in- 
sensitive to possible trends in the data. Simple examples are derivatives of a 
Gaussian, ipQl^gsit) = ^exp(— like the Mexican hat wavelet —"ipcluss 
and the Haar wavelet, ^'nlrl^) = +1 if < t < 1, -1 if 1 < t < 2, and oth- 
erwise. It is straightforward to construct Haar wavelet that are orthogonal to 
linear, quadratic and cubic trends, e. g., "ip^^ls^-it) = 1 for t G [0, 1) U [2, 3), —2 
for t G [1, 2), and otherwise, or i^^Lit) = 1 for t G [0, 1), -3 for t G [1, 2), 
+3 for t G [2,3), -1 for t G [3,4), and otherwise. 



5.2 Discrete Wavelet Transform (WT) Approach 

A detrending fractal analysis of time series can be easily implemented by 



considering Haar wavelet coefficients of the profile Y{j), Eq. (10) [291 [T7]. In 



this case the convolution (16) corresponds to the addition and subtraction 
of mean values of Y{j) within segments of size s. Hence, defining l^(s) = 
s ^j=i yi^s + j), the coefficients can be written as 

FwToii^, s) = L (0) (z/s, s) = Y^{s) - Fj,+i(s), (17) 



laar 



FwTi(i^, s) = L (1) (z/s,s) = Y^{s) - 2Y^+i{s) + F,,+2(s), and (18) 

^Haar 

FwT2(i^, s) = L (2) (z/s, s) = Y,{s) - 3n+i(s) + 3^+2(5) - ^+3(5) (19) 



-laar 



for constant, linear and quadratic detrending, respectively. The generaliza- 
tion for higher orders of detrending is obvious. The resulting mean-square 
fluctuations F^rp^(i/, s) are averaged over all u to obtain the mean fluctua- 



tion F2{s), see Eq. (13). Figure 5 shows typical results for WT analysis of 
long-term correlated, short-term correlated and uncorrelated data. 

Regarding trend-elimination, wavelet transform WTO corresponds to stan- 
dard FA (see Section 4.4), and only constant trends in the proflle are elim- 
inated. WTl is similar to Hurst's rescaled range analysis (see Section 4.3): 
linear trends in the proflle and constant trends in the data are eliminated, 
and the range of the fluctuation exponent a ~ if is up to 2. In general, 
WTn determines the fluctuations from the nth derivative, this way eliminat- 
ing trends described by [n — l)st-order polynomials in the data. The results 
become statistically unreliable for scales s larger than one tenth of the length 
of the data, just as for FA. 



5.3 Detrended Fluctuation Analysis (DFA) 

In the last 13 years Detrended Fluctuation Analysis (DFA), originally intro- 
duced by Peng et al. [30], has been established as an important method to 
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reliably detect long-range (auto-) correlations in non-stationary time series. 
The method is based on random walk theory and basically represents a lin- 
ear detrending version of FA (see Section 4.4). DFA was later generalised for 
higher order detrending [15], separate analysis of sign and magnitude series 
[3T] (see Section 5.5), multifractal analysis |32] (see Section 6.3), and data 
with more than one dimension [33]. Its features have been studied in many 
articles [311 ESI |36l |37l [38l [39]. In addition, several comparisons of DFA 
with other methods for stationary and non-stationary time-series analysis 
have been published, see, e. g., [2T1 SDl SI] and in particular [22], where 
DFA is compared with many other established methods for short data sets, 
and [12], where it is compared with recently suggested improved methods. 
Altogether, there are about 450-500 papers applying DFA (till April 2008). 
In most cases positive auto-correlations were reported leaving only a few 
exceptions with ant i- correlations, see, e. g., [121 SH SS] . 

Like in the FA method, one first calculates the global profile according to 



Eq. (10) and divides the profile into Ng = int(A^/s) non-overlapping segments 
of size s starting from the beginning and another Ns segments starting from 
the end of the considered series. DFA explicitly deals with monotonous 
trends in a detrending procedure. This is done by estimating a polynomial 
trend y^sU) within each segment u by least-square fitting and subtracting 
this trend from the original profile ('detrending'), 

Uj) = Y{j)~yZ{j)- (20) 

The degree of the polynomial can be varied in order to eliminate constant 
(m = 0), linear (m = 1), quadratic (m = 2) or higher order trends of the 
profile function [15] . Conventionally the DFA is named after the order of the 
fitting polynomial (DFAO, DFAl, DFA2, ...). In DFAm, trends of order m in 
the profile Y{j) and of order m — 1 in the original record eliminated. 
The variance of the detrended profile Ys{j) in each segment yields the 
mean-square fiuctuations, 

i^DFA™(^,^) = -En'(j)- (21) 

As for FA and discrete wavelet analysis, the -FoFAml'^) ^) averaged over all 



segments u to obtain the mean fiuctuations -^2(5), see Eq. (14). Calculating 
^2(5) for many s, the fiuctuation scaling exponent a can be determined just 
as with FA. Figure 6 shows typical results for DFA of the same long-term 
correlated, short-term correlated and uncorrelated data studied already in 
Fig. 5. 
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We note that in studies that include averaging over many records (or one 
record cut into many separate pieces by the ehmination of some unrehable 
intermediate data points) the averaging procedure (13) must be performed 
for all data. Taking the square root is always the final step after all averaging 
is finished. It is not appropriate to calculate -^2(5) for parts of the data and 
then average the -^2(5) values, since such a procedure will bias the results 
towards smaller scaling exponents on large time scales. 

If F2{s) increases for increasing s by F2{s) ~ with 0.5 < a < 1, one 
finds that the scaling exponent a ^ H is related to the correlation exponent 7 
by a = 1— 7/2 (see Eq. (15)). A value of a = 0.5 thus indicates that there are 
no (or only short-range) correlations. If a > 0.5 for all scales s, the data are 
long-term correlated. The higher a, the stronger the correlations in the signal 
are. a > 1 indicates a non-stationary local average of the data; in this case, 
FA fails and yields only a = 1. The case a < 0.5 corresponds to long-term 
anti-correlations, meaning that large values are most likely to be followed 
by small values and vice versa, a values below er not possible. Since the 
maximum value for a in DFAm is m + 1, higher detrending orders should 
be used for very non-stationary data with large a. Like in FA and Hurst's 
analysis, a will decrease or increase by one upon additional differentiation or 
integration of the data, respectively. 

Small deviations from the scaling law (14), i. e. deviations from a straight 
line in a double logarithmic plot, occur for small scales s, in particular for 
DFAm with large detrending order m. These deviations are intrinsic to the 
usual DFA method, since the scaling behaviour is only approached asymptot- 
ically. The deviations limit the capability of DFA to determine the correct 
correlation behaviour in very short records and in the regime of small s. 
DFA6, e. g., is only defined for s > 8, and significant deviations from the 
scaling law -^2(5) ~ 5°" occur even up to s ~ 30. They will lead to an over- 
estimation of the fluctuation exponent a, if the regime of small s is used in 
a fitting procedure. An approach for correction of this systematic artefact in 
DFA is described in |34j . 

The number of independent segments of length s is larger in DFA than in 
WT, and the fluctuations in FA are larger than in DFA. Hence, the analysis 
has to be based on s values lower than Smax = N/A for DFA compared 
with Smax = A^/10 for FA and WT. The accuracy of scaling exponents a 
determined by DFA was recently studied as a function of the length of 
the data |12] (fitting range s G [10,A^/2] was used). The results show that 
statistical standard errors of a (one standard deviation) are approximately 
0.1 for = 500, 0.05 for = 3 000, and reach 0.03 for = 10 000. Findings 
of long-term correlations with a = 0.6 in data with only 500 points are thus 
not significant; and a should be at least 0.55 even for data of 10 000 points. 
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A generalization of DFA for two-dimensional data (or even higher dimen- 
sions d) was recently suggested [33]. The generalization works well when 
tested with synthetic surfaces including fractional Brownian surfaces and 
multifractal surfaces. The two-dimensional MFDFA is also adopted to anal- 
yse two images from nature and experiment, and nice scaling laws are unrav- 
elled. In the 2d procedure, a double cumulative sum (profile) is calculated 
by summing over both directional indices analogous with Eq. (10), Y{k, I) = 
Y.i=i I]j=i This surface is partitioned into squares of size s x s with in- 
dices u and fi, in which polynomials like vl^^^sihj) = ai'^+bj'^2+cij+di+ej+f 
are fitted. The fluctuation function -^2(5) is again obtained by calculating 
the variance of the profile from the fits. 



5.4 Detection of Trends and Crossovers with DFA 

Frequently, the correlations of recorded data do not follow the same scaling 
law for all time scales s, but one or sometimes even more crossovers between 
different scaling regimes are observed (see Section 3.3). Time series with a 
well-defined crossover at Sx and vanishing correlations above Sx are most eas- 
ily generated by Fourier filtering (see Section 8.1). The power spectrum S{f) 
of an uncorrelated random series is multiplied by (///x)~^ with j3 = 2a — 1 
for frequencies f > fx = ^/sx only. The series obtained by inverse Fourier 
transform of this modified power spectrum exhibits power-law correlations 
on time scales s < Sx only, while the behaviour becomes uncorrelated on 
larger time scales s > Sx. 

The crossover from F2{s) ~ to -^2(5) ~ s'/^ is clearly visible in double 
logarithmic plots of the DFA fluctuation function for such short-term cor- 
related data. However, it occurs at times Sx*"^ that are different from the 
original Sx used for the generation of the data and that depend on the de- 
trending order m. This systematic deviation is most significant in the DFAm 
with higher m. Extensive numerical simulations (see Fig. 3 in [34j) show that 
the ratios of s^x^V^^x are 1.6, 2.6, 3.6, 4.5, and 5.4 for DFAl, DFA2, . . . , DFA5, 
with an error bar of approximately 0.1. Note, however, that the precise value 
of this ratio will depend on the method used for fitting the crossover times 
Sx"^^ (and the method used for generating the data if generated data is anal- 
ysed). If results for different orders of DFA shall be compared, an observed 
crossover Sx'"'' can be systematically corrected dividing by the ratio for the 
corresponding DFAm. If several orders of DFA are used in the procedure, 
several estimates for the real Sx will be obtained, which can be checked for 
consistency or used for an error approximation. A real crossover can thus 
be well distinguished from the effects of non-stationarities in the data, which 



19 



lead to a different dependence of an apparent crossover on m. 

The procedure is also required if the characteristic time scale of short- 
term correlations shall be studied with DFA. If consistent (corrected) Sx 
values are obtained based on DFAm with different m, the existence of a real 
characteristic correlation time scale is positively confirmed. Note that lower 
detrending orders are advantageous in this case, since the observed crossover 
time scale Sx™^ might become quite large and nearly reach one forth of the 
total series length (iV/4), where the results become statistically inaccurate. 

We would like to note that studies showing scaling long-term correla- 
tions should not be based on DFA or variants of this method alone in most 
applications. In particular, if it is not clear whether a given time series is 
indeed long-term correlated or just short-term correlated with a fairly large 
crossover time scale, results of DFA should be compared with other meth- 
ods. For example, one can employ wavelet methods (see, e. g.. Section 5.2). 
Another option is to remove short-term correlations by considering averaged 
series for comparison. For a time series with daily observations and possible 
short-term correlations up to two years, for example, one might consider the 
series of two-year averages and apply DFA together with FA, binned power 
spectra analysis, and/or wavelet analysis. Only if these methods still indicate 
long-term correlations, one can be sure that the data are indeed long-term 
correlated. 

As discussed in Section 3.3, records from real measurements are often 

affected by non-stationarities, and in particular by trends. They have to be 
well distinguished from the intrinsic fluctuations of the system. To investigate 
the effect of trends on the DFAm fluctuation functions, one can generate 
artificial series {xi) with smooth monotonous trends by adding polynomials 
of different power p to the original record (xi), 

Xi = Xi + AxP with x = i/N. (22) 

For the DFAm, such trends in the data can lead to an artificial crossover 
in the scaling behaviour of i^2(s), i. e., the slope a is strongly increased for 
large time scales s. The position of this artificial crossover depends on the 
strength A and the power p of the trend. Evidently, no artificial crossover 
is observed, if the detrending order m is larger than p and p is integer. The 
order p of the trends in the data can be determined easily by applying the 
different DFAm. If p is larger than m or p is not an integer, an artiflcial 
crossover is observed, the slope atrcnd in the large s regime strongly depends 
on m, and the position of the artiflcial crossover also depends strongly on m. 
The artificial crossover can thus be clearly distinguished from real crossovers 
in the correlation behaviour, which result in identical slopes a and rather 
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similar crossover positions for all detrending orders m. For more extensive 
studies of trends with non-integer powers we refer to [Ml [35] . The effects of 
periodic trends are also studied in [3^ . 

If the functional form of the trend in given data is not known a priori, the 
fluctuation function -^2(3) should be calculated for several orders m of the 
fitting polynomial. If m is too low, F2{s) will show a pronounced crossover 
to a regime with larger slope for large scales s [SU [35] . The maximum slope 
of logF2(s) versus logs is m + 1. The crossover will move to larger scales 
s or disappear when m is increased, unless it is a real crossover not due to 
trends. Hence, one can find m such that detrending is sufficient. However, m 
should not be larger than necessary, because shifts of the observed crossover 
time scales and deviations on short scales s increase with increasing m. 

5.5 Sign and Magnitude (Volatility) DFA 

To study the origin of long-term fractal correlations in a time series, the series 
can be split into two parts, which are analysed separately. It is particularly 
useful to split the series of increments, Axj = Xi — i = 1, . . . , N, into a 
series of signs Xi = Si = sign Ax j and a series of magnitudes Xi = rrii = \Axi\ 
[3TI H6l There is an extensive interest in the magnitude time series 

in economics [HI SHj. These data, usually called volatility, represents the 
absolute variations in stock (or commodity) prices and are used as a measure 
quantifying the risk of investments. While the actual prices are only short- 
term correlated, long-term correlations have been observed in volatility series 

Time series having identical distributions and long-range correlation prop- 
erties can exhibit quite different temporal organizations of the magnitude and 
sign sub-series. The DFA method can be applied independently to both of 
these series. Since in particular the signs are often rather strongly anti- 
correlated and DFA will give incorrect results if a is too close to zero, one 
often studies integrated sign and magnitude series. As mentioned above, 
integration Xi S}=i % increases a by one. 

Most published results report short-term anti-correlations and no long- 
term correlations in the sign series, i. e., Osign < 1/2 for the non-integrated 
signs Si (or agign < 3/2 for the integrated signs) on low time scales and 
«sign — > 1/2 asymptotically for large s. The magnitude series, on the other 
hand, are usually either uncorrelated Omagn = 1/2 (or 3/2) or positively long- 
term correlated ttmagn > 1/2 (or 3/2). It has been suggested that findings 
of ctmagn > 1/2 are related with nonlinear properties of the data and in 
particular multifractality [211 SSI HZ] , if a < 1.5 in standard DFA. Specifically, 
the results suggest that the correlation exponent of the magnitude series 
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is a monotonically increasing function of the multifractal spectrum (i. e., 
the singularity spectrum) width of the original series (see Section 6.1). On 
the other hand, the sign series mainly relates to linear properties of the 
original series. At small time scales s < 16 the standard a is approximately 
the average of Osign and Omagn, if integrated sign and magnitude series are 
analysed. For a > 1.5 in the original series, the integrated magnitude and 
sign series have approximately the same two-point scaling exponents |16] . 
An analytical treatment is presented in |47j . 



5.6 Further Detrending Approaches 

A possible drawback of the DFA method is the occurrence of abrupt jumps 



in the detrended profile Ys{j) (Eq. (20)) at the boundaries between the seg- 



ments, since the fitting polynomials in neighbouring segments are not related. 
A possible way to avoid these jumps would be the calculation of F2{s) based 
on polynomial fits in overlapping windows. However, this is rather time con- 
suming due to the polynomial fit in each segment and is consequently not 
done in most applications. To overcome the problem of jumps several modi- 
fications and extensions of the FA and DFA methods have been suggested in 
the last years. These methods include 

• the Detrended Moving Average technique [301 ED E2] , which we denote 
by Backward Moving Average (BMA) technique (following 



the Centred Moving Average Average (CMA) method [33], an essen- 
tially improved version of BMA, 

the Modified Detrended Fluctuation Analysis (MDFA) [51], which is 
essentially a mixture of old FA and DFA, 

the continuous DFA (CDFA) technique [551 ES], which is particularly 
useful for the detection of crossovers, 



the Fourier DFA [57j . 

a variant of DFA based on empirical mode decomposition (EMD) [58] , 

a variant of DFA based on singular value decomposition (SVD) 
and 

a variant of DFA based on high-pass filtering [61] . 
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Detrended Moving Average techniques will be thoroughly described and 
discussed in the next section. A study comparing DFA with CMA and MDFA 
can be found in |12]. For studies comparing DFA and BMA, see [221 ES]; note 
that [63] also discusses CMA. 

The method we denote as Modified Detrended Fluctuation Analysis (MDFA ) 
|54j . eliminates trends similar to the DFA method. A polynomial is fitted to 
the profile function Y{j) in each segment u and the deviation between the 
profile function and the polynomial fit is calculated, = y{j) — vl^sU) 



(Eq. (20 )). To estimate correlations in the data, this method uses a derivative 



of Fs(j), obtained for each segment z/, by Ay's(j) = Ys{,i + s /2) — Ys{j). Hence, 



the fluctuation function (compare with Eqs. (13) and (21)) is calculated as 
follows: 

1/2 

(23) 



1 AT ^ 

j.T.{ys{3 + s/2)-Uj) 



As in case of DFA, MDFA can easily be generalised to remove higher order 
trends in the data. Since the fitting polynomials in adjacent segments are not 
related, Ys{j) shows abrupt jumps on their boundaries as well. This leads to 
fluctuations of F2{s) for large segment sizes s and limits the maximum usable 
scale to s < N/A as for DFA. The detection of crossovers in the data, however, 
is more exact with MDFA (compared with DFA), since no correction of the 
estimated crossover time scales seems to be needed |12]. 

The Fourier- detrended fluctuation analysis [STj aims to eliminate slow os- 
cillatory trends which are found especially in weather and climate series due 
to seasonal influences. The character of these trends can be rather periodic 
and regular or irregular, and their influence on the detection of long-range 
correlations by means of DFA was systematically studied previously [S^ . 
Among other things it has been shown that slowly varying periodic trends 
disturb the scaling behaviour of the results much stronger than quickly os- 
cillating trends and thus have to be removed prior to the analysis. In case of 
periodic and regular oscillations, e. g., in temperature fluctuations one simply 
removes the low frequency seasonal trend by subtracting the daily mean tem- 
peratures from the data. Another way, which the Fourier-detrended fluctua- 
tion analysis suggests, is to filter out the relevant frequencies in the signals' 
Fourier spectrum before applying DFA to the filtered signal. Nevertheless, 
this method faces several difficulties especially its limitation to periodic and 
regular trends and the need for a priori knowledge of the interfering frequency 
band. 

To study correlations in data with quasi-periodic or irregular oscillating 
trends, empirical mode decomposition (EMD) was suggested [58|- The EMD 
algorithm breaks down the signal into its intrinsic mode functions (IMFs) 
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which can be used to distinguish between fluctuations and background. The 
background, estimated by a quasi-periodic fit containing the dominating fre- 
quencies of a sufficiently large number of IMFs, is subtracted from the data, 
yielding a slightly better scaling behaviour in the DFA curves. However, we 
believe that the method might be too complicated for wide-spread applica- 
tions. 

Another method which was shown to minimise the effect of periodic and 
quasi-periodic trends is based on singular value decomposition (SVD) [H^ISU]. 
In this approach, one first embeds the original signal in a matrix whose di- 
mension has to be much larger than the number of frequency components of 
the periodic or quasi-periodic trends obtained in the power spectrum. Ap- 
plying SVD yields a diagonal matrix which can be manipulated by setting 
the dominant eigen-values (associated with the trends) to zero. The filtered 
matrix finally leads to the filtered data, and it has been shown that subse- 
quent application of DFA determines the expected scaling behaviour if the 
embedding dimension is sufficiently large. None the less, the performance of 
this rather complex method seems to decrease for larger values of the scaling 
exponent. Furthermore SVD-DFA assumes that trends are deterministic and 
narrow banded. 

The detrending procedure in DFA (Eq. (20)) can be regarded as a scale- 
dependent high-pass filter since (low-frequency) fluctuations exceeding a spe- 
cific scale s are eliminated. Therefore, it has been suggested to obtain the 
detrended profile Ys{j) for each scale s directly by applying digital high- 
pass filters |61j- In particular, Butterworth, Chebyshev-I, Chebyshev-II, 
and an elliptical filter were suggested. While the elliptical filter showed the 
best performance in detecting long-range correlations in artificial data, the 
Chebyshev-II filter was found to be problematic. Additionally, in order to 
avoid a time shift between filtered and original profile, the average of the di- 
rectly filtered signal and the time reversed filtered signal is considered. The 
effects of these complicated filters on the scaling behaviour are, however, not 
fully understood. 

Finally, a continuous DFA method has been suggested in the context 
of studying heartbeat data during sleep [55t 156] . The method compares 
unnormalised fluctuation functions F2{s) for increasing length of the data. 
I. e., one starts with a very short recording and subsequently adds more 
points of data. The method is particularly suitable for the detection of 
change points in the data, e. g., physiological transitions between different 
activity or sleep stages. Since the main objective of the method is not the 
study of scaling behaviour, we do not discuss it in detail here. 
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5.7 Centered Moving Average (CMA) Analysis 

Particular attractive modifications of DFA are tlie Detrended Moving Average 
(DMA) metfiods, wfiere running averages replace tlie polynomial fits. Tlie 
first suggested version, the Backward Moving Average (BMA) method [50l 
|511[52], however, suffers from severe problems, because an artificial time shift 
of s between the original signal and the moving average is introduced. This 
time shift leads to an additional contribution to the detrended profile 
which causes a larger fluctuation function ^2(5) in particular for small scales 
in the case of long-term correlated data. Hence, the scaling exponent a is 
systematically underestimated [62]. In addition, the BMA method preforms 
even worse for data with trends [63], and its slope is limited by a < 1 just 
as for the non-detrending method FA. 

It was soon recognised that the intrinsic error of BMA can be overcome 
by eliminating the artificial time shift. This leads to the Centred Moving 
Average (CMA) method [53], where Ys{j) is calculated as 

Ysij)=Yij)-- ^'2^' F(j+z), (24) 

^ i=-(s-l)/2 



replacing Eq. (20) while Eq. (21 ) and the rest of the DFA procedure described 
in Section 5.3 stay the same. Unlike DFA, the CMA method cannot easily 
be generalised to remove linear and higher order trends in the data. 

It was recently proposed [12] that the scaling behaviour of the CMA 
method is more stable than for DFAl and MDFAl, suggesting that CMA 
could be used for reliable computation of a even for scales s < 10 (without 
correction of any systematic deviations needed in DFA for this regime) and 
up to Smax = A^/2. The standard errors in determining the scaling exponent 
a by fitting straight lines to the double logarithmic plots of F2{s) have been 
studied in [12]; they are comparable with DFAl (see end of Section 5.3). 

Regarding the determination of crossovers, CMA is comparable to DFAl. 
Ultimately, the CMA seems to be a good alternative to DFAl when analysing 
the scaling properties in short data sets without trends. Nevertheless for data 
with possible unknown trends we recommend the application of standard 
DFA with several different detrending polynomial orders in order to distin- 
guish real crossovers from artificial crossovers due to trends. In addition, 
an independent approach (e. g., wavelet analysis) should be used to confirm 
findings of long-term correlations (see also Section 5.4). 
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6 Methods for Mult ifr act al Time Series Anal- 
ysis 

This chapter describes the multifractal characterisation of time series, for an 
introduction, see Section 3.4. The simplest type of multifractal analysis is 
based upon the standard partition function multifractal formalism, which has 
been developed for the multifractal characterisation of normalised, stationary 
measures [HI [121 EH [65] . Unfortunately, this standard formalism does not give 
correct results for non-stationary time series that are affected by trends or 
that cannot be normalised. Thus, in the early 1990s an improved multifrac- 
tal formalism has been developed, the wavelet transform modulus maxima 
(WTMM) method [Ml [HZl IMl ISHl [ZO], which is based on wavelet analysis 
and involves tracing the maxima lines in the continuous wavelet transform 
over all scales. An important alternative is the multifractal DFA (MF-DFA) 
algorithm [32], which does not require the modulus maxima procedure, and 
hence involves little more effort in programming than the conventional DFA. 
For studies comparing methods for detrending multifractal analysis (multi- 
fractal DFA (MF-DFA) and wavelet transform modulus maxima (WTMM) 
method), see [32l [711 [12] ■ 

6.1 The Structure Function Approach and Singularity 
Spectra 

In the general multifractal formalism, one considers a normalised measure 
/x(t), t e [0,1], and defines the box probabilities jls{t) = Jl:l^^l2 fi{t') dt' in 
neighbourhoods of (scale) length s -C 1 around t. The multifractal approach 
is then introduced by the partition function 

l/s-l 

Z,(s) = ^ /i^[(z/ + l/2)s] ~s"(<?) for s<l, (25) 

where r(g) is the Renyi scaling exponent and g is a real parameter that 
can take positive as well as negative moments. Note that r(g) is sometimes 
defined with opposite sign (see, e. g., [H]). A record is called monofractal 
(or self-affine), when the Renyi scaling exponent r(g) depends linearly on g; 
otherwise it is called multifractal. The generalised multifractal dimensions 
D{q) (see also Section 3.4) are related to r(g) by D{q) = T{q)/q— 1, such 
that the fractal dimension of the support is -D(O) = — t(0) and the correlation 
dimension is D{2) = t(2). 

In time series, a discrete version has to be used, and the considered data 
(xi), i = 1, . . . ,N may usually include negative values. Hence, setting Ng = 
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int(A^/ s) and s) = J2i=i Xus+i ioi u = 0, . . . , Ng — 1 we can define [6l [12] , 



Ns-l 
:/=0 



for s > 1. 



(26) 



Inserting tfie profile Y{j) and -Ffa(z^, s) from Eqs. (10) and (11), respectively, 
we obtain 



N,-l 



Z,{s)= E {[Y{{u + l)s)-Y{us] 



v=0 



Ns-l 



(27) 



Comparing Eqs. (27) with (13), we see that this multifractal approach can be 
considered as a generalised version of the Fluctuation Analysis (FA) method, 
where the exponent 2 is replaced by q. In particular we find (disregarding 
the summation over the second partition of the time series) 



Fois] 



1/2 



~ s 



[l+r(2)]/2 



2a = 1 + r(2) = 1 + D(2). (28) 



We thus see that all methods for (mono-)fractal time analysis (discussed in 
Chapters 4 and 5) in fact study the correlation dimension D{2) = 2a — 1 = 
/3 = 1 - 7 (see Eq. 

It is straightforward to define a generalised (multifractal) Hurst exponent 
h{q) for the scaling behaviour of the qth moments of the fluctuations [Mf 165]. 



■Zois) 



l/<7 



~ S 



[i+T{g)]/g 



Mq) 



h{q) 



1 + r{q) 



(29) 



with h{2) = a ^ H. In the following, we will use only h{2) for the standard 
fluctuation exponent (denoted by a in the previous chapters), and reserve 
the letter a for the Holder exponent. 

Another way to characterise a multifractal series is the singularity spec- 
trum /(a), that is related to r(g) via a Legendre transform [6l [12], 



a 



d 
dq 



T{q) and /(a) = qa — T{q). 



(30) 



Here, a is the singularity strength or Holder exponent (see also articles ... 
in the encyclopedia), while /(a) denotes the dimension of the subset of the 
series that is characterised by a. Note that a is not the fluctuation scaling 
exponent in this section, although the same letter is traditionally used for 
both. Using Eq. ([29]), we can directly relate a and /(a) to h{q), 



a = h{q) + qh\q) and /(a) = q[a — h{q)] + 1. 



(31) 
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6.2 Wavelet Transform Modulus Maxima (WTMM) 
Method 

The wavelet transform modulus maxima (WTMM) method [66 | [67ll68t |69 | [70] 
is a well-known method to investigate the multifractal scaling properties 
of fractal and self-affine objects in the presence of non-stationarities. For 
applications, see e. g. [73], |71]. It is based upon the wavelet transform with 



continuous basis functions as defined in Section 5.1, Eq. (16). Note that 



in this case the series Xi are analysed directly instead of the profile Y{j) 



defined in Eq. (10). Using wavelets orthogonal to mth order polynomials, 
the corresponding trends are elliminated. 

Instead of averaging over all wavelet coefficients L^{t,s), one averages, 
within the modulo-maxima method, only the local maxima of |L^(r, s)|. 
First, one determines for a given scale s, the positions tj of the local max- 
ima of |Vr(r, s)| as function of r, so that \L^{Tj — l,s)| < \L^{Tj,s)\ > 
\L^{Tj + 1, s)| for j = 1, . . . , jmax- This maxima procedure is demonstrated 
in Fig. 7. Then one sums up the gth power of the maxima, 

Jmax 

Z{q,s) = Y.\L4r,,s)\''. (32) 



The reason for the maxima procedure is that the absolute wavelet coefficients 
|L^(r, s)| can become arbitrarily small. The analysing wavelet ipi^) must 
always have positive values for some x and negative values for other x, since 
it has to be orthogonal to possible constant trends. Hence there are always 



positive and negative terms in the sum (16), and these terms might cancel. 
If that happens, \L^{t, s)\ can become close to zero. Since such small terms 



would spoil the calculation of negative moments in Eq. (32), they have to be 
eliminated by the maxima procedure. 

In fluctuation analysis, on the other hand, the calculation of the variances 
F^(z/, s), e. g. in Eq. (11), involves only positive terms under the summa- 



tion. The variances cannot become arbitrarily small, and hence no maximum 
procedure is required for series with compact support. In addition, the vari- 
ances will always increase if the segment length s is increased, because the 
fit will always be worse for a longer segment. In the WTMM method, in 
contrast, the absolute wavelet coefficients |L^(r, s)| need not increase with 
increasing scale s, even if only the local maxima are considered. The val- 
ues \L^{t, s) \ might become smaller for increasing s since just more (positive 



and negative) terms are included in the summation (16), and these might 
cancel even better. Thus, an additional supremum procedure has been in- 
troduced in the WTMM method in order to keep the dependence of Z{q, s) 
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on s monotonous. If, for a given scale s, a maximum at a certain position tj 
happens to be smaller than a maximum at rj ^ Tj for a lower scale s' < s, 
then L^{Tj, s) is replaced by L^{tj, s') in Eq. (32). 



Often, scaling behaviour is observed for Z(q,s), and scaling exponents 
f{q) can be defined that describe how Z{q,s) scales with s, 

Z{q,s) s^'~'^\ (33) 

The exponents f{q) characterise the multifractal properties of the series under 
investigation, and theoretically they are identical with the r(g) defined in 



Eq. (26) [661 EH EHl CD] and related to h{q) by Eq. (29) 



6.3 Multifractal Detrended Fluctuation Analysis (MF- 
DFA) 

The multifractal DFA (MF-DFA) procedure consists of five steps [32]. The 
first three steps are essentially identical to the conventional DFA procedure 
(see Section 5.3 and Fig. 4). Let us assume that (xj) is a series of length N, 
and that this series is of compact support. The support can be defined as 
the set of the indices j with nonzero values Xj , and it is compact if xj = for 
an insignificant fraction of the series only. The value of xj = is interpreted 
as having no value at this j. Note that we are not discussing the fractal 
or multifractal features of the plot of the time series in a two-dimensional 
graph (see also the discussion in Section 3.1), but analysing time series as 
one- dimensional structures with values assigned to each point. Since real 
time series always have finite length A^, we explicitly want to determine the 
multifractality of finite series, and we are not discussing the limit for N ^ oo 
here (see also Section 6.1). 



Step 1: Calculate the profile Y{j, Eq. (10), by integrating the time 
series. 

Step 2: Divide the profile Y{j) into A^^ = int(A^/s) non-overlapping 
segments of equal length s. Since the length A^ of the series is often not 
a multiple of the considered time scale s, the same procedure can be 
repeated starting from the opposite end. Thereby, 2A^5 segments are 
obtained altogether. 

Step 3: Calculate the local trend for each of the 2Ns segments by a 



least-square fit of the profile. Then determine the variance by Eqs. (20 ) 
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and (21) for each segment z/ = 0, . . . , 2Ns — 1 Again, linear, quadratic, 
cubic, or higher order polynomials can be used in the fitting procedure, 
and the corresponding methods are thus called MF-DFAl, MF-DFA2, 
MF-DFA3, . . .) [32]. In (MF-)DFAm [mth order (MF-)DFA] trends of 
order mp in the profile (or, equivalently, of order m — 1 in the original 
series) are eliminated. Thus a comparison of the results for different 
orders of DFA allows one to estimate the type of the polynomial trend 
in the time series 



Step 4'- Average over all segments to obtain the gth order fluctuation 
function 

f 1 r 10/2^^^ 

F,{s) = \^^^j:[FUmi'^,^)] j , (34) 



This is the generalization of Eq. ( 13 ) suggested by the relations derived 



in Section 6.1. For q = 2, the standard DFA procedure is retrieved. 
One is interested in how the generahsed q dependent fluctuation func- 
tions Fq{s) depend on the time scale s for different values of q. Hence, 
we must repeat steps 2 to 4 for several time scales s. It is apparent 
that Fg{s) will increase with increasing s. Of course, Fq{s) depends on 
the order m. By construction, Fg{s) is only defined for s > m + 2. 



• Step 5: Determine the scaling behaviour of the fluctuation functions by 
analysing log- log plots Fq{s) versus s for each value of q. If the series 
Xi are long-range power-law correlated, Fq[s) increases, for large values 
of s, as a power-law, 

Fq{s) ~ with h{q) = l±lM. (35) 



For very large scales, s > N/A, Fq{s) becomes statistically unreliable 
because the number of segments Ns for the averaging procedure in step 4 
becomes very small. Thus, scales s > N/4 should be excluded from the 
fitting procedure determining h{q). Besides that, systematic deviations from 
the scaling behaviour in Eq. (35), which can be corrected, occur for small 
scales s ~ 10. 

The value of h{0), which corresponds to the limit h{q) for q —>■ 0, cannot 
be determined directly using the averaging procedure in Eq. (34) because of 
the diverging exponent. Instead, a logarithmic averaging procedure has to 
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be employed, 

Fo{s) = exp | — E In [f\u, s)] | ~ .'^W. (36) 

Note that h{0) cannot be defined for time series with fractal support, where 
h{q) diverges for g — 0. 

For monofractal time series with compact support, h{q) is independent 
of q, since the scaling behaviour of the variances F^pAmi^y ^) identical for 
all segments z/, and the averaging procedure in Eq. (34) will give just this 
identical scaling behaviour for all values of q. Only if small and large fluctu- 
ations scale differently, there will be a significant dependence of h{q) on q: If 
we consider positive values of q, the segments u with large variance -F^(z/, s) 
(i. e. large deviations from the corresponding fit) will dominate the average 
Fq{s). Thus, for positive values of q, h{q) describes the scaling behaviour of 
the segments with large fluctuations. On the contrary, for negative values of 
q, the segments u with small variance -FpFAml'^) ^) ^ill dominate the average 
Fq{s). Hence, for negative values of q, h{q) describes the scaling behaviour of 
the segments with small fluctuations. Figure 8 shows typical results obtained 
for Fg{s) in the MF-DFA procedure. 

Usually the large fluctuations are characterised by a smaller scaling ex- 
ponent h{q) for multifractal series than the small fluctuations. This can be 
understood from the following arguments: For the maximum scale s = N the 
fluctuation function Fq{s) is independent of q, since the sum in Eq. (34) runs 
over only two identical segments. For smaller scales s <^ N the averaging 
procedure runs over several segments, and the average value Fq{s) will be 
dominated by the F'^{i', s) from the segments with small (large) fluctuations 
if g < (g > 0). Thus, for s N, Fq{s) with g < will be smaller than 
Fq{s) with g > 0, while both become equal for s = N. Hence, if we assume 
an homogeneous scaling behaviour of Fq{s) following Eq. (35), the slope h{q) 
in a log- log plot of Fq{s) with g < versus s must be larger than the cor- 
responding slope for Fq{s) with g > 0. Thus, h{q) for g < will usually be 
larger than h{q) for g > 0. 

However, the MF-DFA method can only determine positive generalised 
Hurst exponents h{q), and it already becomes inaccurate for strongly anti- 
correlated signals when h{q) is close to zero. In such cases, a modified (MF- 
)DFA technique has to be used. The most simple way to analyse such data 
is to integrate the time series before the MF-DFA procedure. Following the 
MF-DFA procedure as described above, we obtain a generalised fluctuation 
functions described by a scaling law with h{q) = h{q) + 1. The scaling 
behaviour can thus be accurately determined even for h{q) which are smaller 
than zero for some values of g. 
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The accuracy of h{q) determined by MF-DFA certainly depends on the 
length of the data. For q = ±10 and data with = 10 000 and 100 000, 
systematic and statistical error bars (standard deviations) up to Ah{q) ^ 
±0.1 and ~ ±0.05 should be expected, respectively [32]. A difference of 
h{—10) — h{+10) = 0.2, corresponding to an even larger width Aa of the 
singularity spectrum /(a) defined in Eq. (30) is thus not significant unless 
the data was longer than A^ = 10 000 points. Hence, one has to be very 
careful when concluding multifractal properties from differences in h{q). 

As already mentioned in the introduction, two types of multifractality 
in time series can be distinguished. Both of them require a multitude of 
scaling exponents for small and large fluctuations: (i) Multifractality of a 
time series can be due to a broad probability density function for the values 
of the time series, and (ii) multifractality can also be due to different long- 
range correlations for small and large fluctuations. The most easy way to 
distinguish between these two types is by analysing also the corresponding 
randomly shuffled series [3^. In the shuffling procedure the values are put 
into random order, and thus all correlations are destroyed. Hence the shuffled 
series from multifractals of type (ii) will exhibit simple random behaviour, 
^shuf(?) = 0.5, i. e. non-multifractal scaling. For multifractals of type (i), on 
the contrary, the original h{q) dependence is not changed, h{q) = hshuiig), 
since the multifractality is due to the probability density, which is not affected 
by the shuffling procedure. If both kinds of multifractality are present in a 
given series, the shuffled series will show weaker multifractality than the 
original one. 



6.4 Comparison of WTMM and MF-DFA 

The MF-DFA results turn out to be slightly more reliable than the WTMM 
results Ell 112] • In particular, the MF-DFA has slight advantages for 
negative q values and short series. In the other cases the results of the two 
methods are rather equivalent. Besides that, the main advantage of the MF- 
DFA method compared with the WTMM method lies in the simplicity of 
the MF-DFA method. However, contrary to WTMM, MF-DFA is restricted 
to studies of data with full one-dimensional support, while WTMM is not. 
Both, WTMM and MF-DFA have been generalised for higher dimensional 
data, see p3] for higher dimensional MF-DFA and, e. g., [70] for higher 
dimensional WTMM. Studies of other generalisations of detrending methods 
like the discrete WT approach (see Section 5.2) and the CMA method (see 
Section 5.7) are currently under investigation [75] . 
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7 Statistics of Extreme Events in Fractal Time 
Series 



The statistics of return intervals between well defined extremal events is a 
powerful tool to characterise the temporal scaling properties of observed time 
series and to derive quantities for the estimation of the risk for hazardous 
events like floods, very high temperatures, or earthquakes. It was shown 
recently that long-term correlations represent a natural mechanism for the 
clustering of the hazardous events [76]. In this chapter we will discuss the 
most important consequences of long-term correlations and fractal scaling of 
time series upon the statistics of extreme events |76l [771 EHl EHl EO]- Corre- 
sponding work regarding multifractal data [HI] is not discussed here. 



7.1 Return Intervals Between Extreme Events 

To study the statistics of return intervals we consider again a time series 
(xi), i = 1,...,N with fractal scaling behaviour, sampled homogeneously 
and normalised to zero mean and unit variance. For describing the reoc- 
currence of rare events exceeding a certain threshold q, we investigate the 
return intervals Vg between these events, see Fig. 9. The average return 
interval Rg = (r^) increases as a function of the threshold q (see, e. g. 
|82j). It is known that for uncorrelated records ('white noise'), the return 
intervals are also uncorrelated and distributed according to the Poisson dis- 
tribution, Pg{r) = exp(— r/i?g). For fractal (long-term correlated) data 
with auto-correlations following Eq. ([5]), we obtain a stretched exponential 
[761I771I781I791E3], 

P,ir) = ^exp[-b,{r/Rgn (37) 

Kg 

This behaviour is shown in Fig. 10. The exponent 7 is the correlation 
exponent from C{s), and the parameters and are independent of q. 
They can be determined from the normalization conditions for Pg{r), i. e.. 



/ Pq{r) dr = 1 and / rPg{r)dr = Rg. The form of the distribution (37) 
indicates that return intervals both well below and well above their average 
value Rg (which is independent of 7) are considerably more frequent for 
long-term correlated than for uncorrelated data. It has to be noted that 
there are deviations from the stretched exponential law ( j37) ) for very small 
r (discretization effects and an additional power-law regime) and for very 
large r (finite size effects), see Fig. 10. The extent of the deviations from 



Eq. (37) depends on the distribution of the values Xi of the time series. For 



a discussion of these effects, see [79] . 
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Equation (37) does not quantify, however, if the return intervals them- 
selves are arranged in a correlated or in an uncorrelated fashion, and if clus- 
tering of rare events may be induced by long-term correlations. To study 
this question, one has to evaluate the auto-correlation function Cr{s) = 
{rqiVjTqil + s)) — B?q oi the return intervals. The results for model data 
suggests that also the return intervals are long-term power-law correlated, 
with the same exponent 7 as the original record. Accordingly, large and 
small return intervals are not arranged in a random fashion but are expected 
to form clusters. As a consequence, the probability of finding a certain return 
interval r depends on the value of the preceding interval ro, and this effect 
has to be taken into account in predictions and risk estimations [761 ES] ■ 

The conditional distribution function Pg(r|ro) is a basic quantity, from 
which the relevant quantities in risk estimations can be derived [82J. For 
example, the first moment of Pg(r|ro) is the average value Rq{ro) of those 
return intervals that directly follow tq. By definition, -Rg(ro) is the expected 
waiting time to the next event, when the two events before were separated 
by tq. The more general quantity is the expected waiting time Tq{x\ro) to 
the next event, when the time x has been elapsed. For x = 0, rg(0|ro) is 
identical to Rqiro). In general, rg(x|ro) is related to -Pg(r|ro) by 

Tq{x\rQ) = I {r — x)Pq{r\ro)dr / / Pq{r\ro)dr (38) 

Jx J X 

For uncorrelated records, Tq{x\ro) / Rq = 1 (except for discreteness effects 
that lead to Tq{x\ro)/Rq > 1 for a; > 0, see [HI])- Due to the scaling of 
Pg(r|ro), also Tq{x\rQ) / Rq scales with ro/Rg and x/ Rq. Small and large return 
intervals are more likely to be followed by small and large ones, respectively, 
and hence rq(0|ro)/-Rg = Rq{ro)/Rq is well below (above) one for r^/Rq well 
below (above) one. With increasing x, the expected residual time to the next 
event increases. Note that only for an infinite long-term correlated record, 
the value of Tq(x|ro) will increase indefinitely with x and tq. For real (finite) 
records, there exists a maximum return interval which limits the values of x, 
ro and rg(x|ro). 



7.2 Distribution of Extreme Events 

In this section we describe how the presence of fractal long-term correlations 
affects the statistics of the extreme events, i. e., maxima within time seg- 
ments of fixed duration i?, see Fig. 11 for illustration. By definition, extreme 
events are rare occurrences of extraordinary nature, such as, e. g. floods, 
very high temperatures, or earthquakes. In hydrological engineering such 
conventional extreme value statistics are commonly applied to decide what 
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building projects are required to protect riverside gainst typical floods 

that occur for example once in 100 years. Most of these results are based 
on statistically independent values Xi and hold only in the limit R oo. 
However, both of these assumptions are not strictly fulfilled for correlated 
fractal scaling data. 

In classical extreme value statistics one assumes that records (xj) consist 
of i.i.d. data, described by density distributions P{x), which can be, e. g., 
a Gaussian or an exponential distribution. One is interested in the distri- 
bution density function Pfj(m) of the maxima {rrij) determined in segments 
of length R in the original series (xj), see Fig. 11. Note that all maxima 
are also elements of the original data. The corresponding integrated maxima 
distribution Gr^tji) is defined as 



Since Gpiim) is the probability of finding a maximum smaller than m, Epi{m) 
denotes the probability of finding a maximum that exceeds m. One of the 
main results of traditional extreme value statistics states that for indepen- 
dently and identically distributed (i.i.d.) data (xj) with Gaussian or expo- 
nential distribution density function P{x) the integrated distribution Gpi{m) 
converges to a double exponential (Fisher-Tippet-Gumbel) distribution (of- 
ten labelled as Type I) (HSl EHl EH EHl EH], i- e.. 



for i? — > oo, where a is the scale parameter and u the location parameter. 
By the method of moments those parameters are given by a = —aji and 
u = mR- n^a with the Euler constant = 0.577216 [881 EQl EH ES] . Here 
rrifi and aji denote the (i?-dependent) mean maximum and the standard 
deviation, respectively. Note that different asymptotics will be reached for 
broader distributions of data (xj) that belong to other domains of attraction 
[88j . For example, for data following a power-law distribution (or Pareto 
distribution), P{x) = {x/xq)~'', Gnim) converges to a Frechet distribution, 
often labelled as Type II. For data following a distribution with finite upper 
endpoint, for example the uniform distribution P{x) = 1 for < a; < 1 , 
Gnirn) converges to a Weibull distribution, often labelled as Type III. We 
do not consider the latter two types of asymptotics here. 

Numerical studies of fractal model data have recently shown that the 
distribution P{x) of the original data has a much stronger effect upon the 
convergence towards the Gumbel distribution than the long-term correla- 
tions in the data. Long-term correlations just slightly delay the convergence 




(39) 




(40) 
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of Gjilm) towards the Gumbel distribution (40). This can be observed very 
clearly in a plot of the integrated and scaled distribution G/j(m) on logarith- 
mic scale [SU] . 

Furthermore, it was found numerically that (i) the maxima series {frij) 
exhibit long-term correlations similar to those of the original data (xj), and 
most notably (ii) the maxima distribution as well as the mean maxima sig- 
nificantly depend on the history, in particular on the previous maximum [80] . 
The last item implies that conditional mean maxima and conditional maxima 
distributions should be considered for improved extreme event predictions. 



8 Simple Models for Fractal and Multifractal 
Time Series 

8.1 Fourier Filtering 

Fractal scaling with long-term correlations can be introduced most easily 
into time series by the Fourier-filtering technique, see, e. g., [931 IHU |95] . The 
Fourier filtering technique is not limited to the generation of long-term corre- 
lated data characterised by a power-law auto-correlation function C{s) ~ x~'^ 
with < 7 < 1. All values of the scaling exponents a = h{2) ^ H or 
P = 2a — 1 can be obtained, even those that cannot be found directly by 
the fractal analysis techniques described in Chapters 4 and 5 (e. g. a < 1). 
Note, however, that Fourier filtering will always yield Gaussian distributed 
data values and that no nonlinear or multifractal properties can be achieved 
(see also Sections 3.4, 5.5, and Chapter 6). In Section 5.4, we have briefly 
described a modification of Fourier filtering for obtaining reliable short-term 
correlated data. 

For the generation of data characterised by fractal scaling with /5 = 2a; — 1 
[93l [9l] we start with uncorrelated Gaussian distributed random numbers 
Xi from an i.i.d. generator. Transforming a series of such numbers into 
frequency space with discrete Fourier transform or FFT (for suitable series 
lengths A^) yields a fiat power spectrum, since random numbers correspond to 
white noise. Multiplying the (complex) Fourier coefficients by /~^''^, where 
/ oc 1/s is the frequency, will rescale the power spectrum S{f) to follow 
Eq. ([6]), as expected for time series with fractal scaling. After transforming 
back to the time domain (using inverse Fourier transform or inverse FFT) 
we will thus obtain the desired long-term correlated data x,. The final step 
is the normalization of this data. 

The Fourier filtering method can be improved using modified Bessel func- 
tions instead of the simple factors f~^^'^ in modifying the Fourier coefficients 
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|95j . This way problems with the divergence of the autocorrelation function 
C(s) at s = can be avoided. 

An alternative method to the Fourier filtering technique, the random mid- 
point displacement method, is based on the construction of self-affine surfaces 
by an iterative procedure, see, e. g. [6]. Starting with one interval with con- 
stant values, the intervals are iterative split in the middle and the midpoint is 
displaced by a random offset. The amplitude of this offset is scaled according 
to the length of the interval. Since the method generates a self-affine surface 
Xi characterised by a Hurst exponent H, the differentiated series Axi can 
be used as long-term correlated or anti-correlated random numbers. Note, 
however, that the correlations do not persist for the whole length of the data 
generated this way. Another option is the use of wavelet synthesis, the reverse 
of wavelet analysis described in Section 5.1. In that method, the scaling law 
is introduced by setting the magnitudes of the wavelet coefficients according 
to the corresponding time scale s. 



8.2 The Schmitz-Schreiber Method 

When long-term correlations in random numbers are introduced by the Fourier- 
filtering technique (see previous section), the original distribution P{x) of the 
time series values Xi is always modified such that it becomes closer to a Gaus- 
sian. Hence, no series (xj) with broad distributions of the values and fractal 
scaling can be generated. In these cases an iterative algorithm introduced by 
Schreiber and Schmitz HE] must be applied. 

The algorithm consists of the following steps: First one creates a Gaus- 
sian distributed long-term correlated data set with the desired correlation 
exponent 7 by standard Fourier-filtering [95]. The power spectrum Scif) = 
-^g(/)-Fg(/) of this data set is considered as reference spectrum (where / de- 
notes the frequency in Fourier space and the -Fg(/) are the complex Fourier 
coefficients). Next one creates an uncorrelated sequence of random numbers 
(xf^), following a desired distribution P{x). The (complex) Fourier trans- 
form F{f) of the (xf^) is now divided by its absolute value and multiplied 
by the square root of the reference spectrum, 

Fif)JSG{f) 

FneM) = \F{f)\ ■ ^^^^ 

After the Fourier back-transformation of F^^ewif), the new sequence {x^^"") 
has the desired correlations (i. e. the desired 7), but the shape of the distri- 
bution has changed towards a (more or less) Gaussian distribution. In order 
to enforce the desired distribution, we exchange the (xf™) by the {xY'^), such 
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that the largest value of the new set is replaced by the largest value of the 
reference set, the second largest of the new set by the second largest of the 
reference set and so on. After this the new sequence has the desired distri- 
bution and is clearly correlated. However, due to the exchange algorithm 
the perfect long-term correlations of the new data sequence were slightly 
altered again. So the procedure is repeated: the new sequence is Fourier 
transformed followed by spectrum adjustment, and the exchange algorithm 
is applied to the Fourier back-transformed data set. These steps are repeated 
several times, until the desired quality (or the best possible quality) of the 
spectrum of the new data series is achieved. 



8.3 The Extended Binomial Multifractal Model 

The multifractal cascade model [6l|6ll|32] is a standard model for multifractal 
data, which is often applied, e. g., in hydrology [96j. In the model, a record Xi 
of length = 2"™="" is constructed recursively as follows. In generation n = 0, 
the record elements are constant, i. e. Xj = 1 for alH = 1, . . . , A^. In the first 
step of the cascade (generation n = 1), the first half of the series is multiplied 
by a factor a and the second half of the series is multiplied by a factor b. 
This yields Xi = a for i = 1, . . . , N/2 and Xi = b for i = N/2 + 1, . . . , N . 
The parameters a and b are between zero and one, < a < 6 < 1. One 
need not restrict the model to 6 = 1 — a as is often done in the literature 
[B]. In the second step (generation n = 2), we apply the process of step 
1 to the two subseries, yielding Xi = a?' for z = 1, . . . , A^/4, Xi = ab for i = 
N/A+1, . . . , A^/2, Xi = ba = ab for i = N/2 + 1, . . . , 3N/A, and Xi = 6^ for i = 
3N/4 + 1, . . . ,N. In general, in step n + 1, each subseries of step n is divided 
into two subseries of equal length, and the first half of the Xi is multiplied by a 
while the second half is multiplied by b. For example, in generation n = 3 the 
values in the eight subseries are a^, a^b, a^b, ab"^, a^b, ab^, ab"^, b^. After 
^max steps, the final generation has been reached, where all subseries have 
length 1 and no more splitting is possible. We note that the final record can 
be written as Xj = a"max-n(«-i)^n(«-i)^ where n{i) is the number of digits 1 in 
the binary representation of the index i, e. g. 72(13) = 3, since 13 corresponds 
to binary 1101. 

For this multiplicative cascade model, the formula for r(g) has been de- 
rived earher P Ell |32] . The result is r(g) = [- ln(a« + 6«) + g ln(a + b)]/ In 2 
or 

^ 1 _ ln(a» + t.) Mo + fO 

q gln2 ln2 ^ ^ 

It is easy to see that h{l) = 1 for all values of a and b. Thus, in this form the 
model is limited to cases where h{l), which is the exponent Hurst defined 
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originally in the R/S method, is equal to one. 

In order to generalise this multifractal cascade process such that any value 
of h{l) is possible, one can subtract the offset A/i = ln(a + 6)/ln(2) from h{q) 
|99j . The constant offset Ah corresponds to additional long-term correlations 
incorporated in the multiplicative cascade model. For generating records 
without this offset, we rescale the power spectrum. First, we fast-Fourier 
transform (FFT) the simple multiplicative cascade data into the frequency 
domain. Then, we multiply all Fourier coefficients by /~^'*, where / is the 
frequency. This way, the slope (3 of the power spectra S{f) ~ is decreased 
from/5 = 2/i(2)-l = [2 ln(a+6)-ln(a2+62)]/ in 2 into /5' = 2[/i(2)-A/i]-l = 
— ln(a^+6^) / In 2. Finally, backward FFT is employed to transform the signal 
back into the time domain. 



8.4 The Bi-fractal Model 

In some cases a simple bi-fractal model is already sufficient for modelling 
apparently multifractal data |100j . For bi-fractal records the Renyi exponents 
r(g) are characterised by two distinct slopes ai and 02, 

-(g) = j 7 - '^ , (43) 

[ ga2 + gx («! - "2) - 1 g > gx 

or 

t{q) = / + ("2 - «i) - 1 g < gx /^^x 



If this behaviour is translated into the h{q) picture using Eq. (29), we obtain 
that h{q) exhibits a plateau from q = —00 up to a certain q^ and decays 
hyperbolically for q > qx, 

%) = |"\'-^i, , , (45) 

I gx(«i - "2)^ + "2 g > gx ^ ^ 



or vice versa, 



= | gx(«2-a0j + «i g<gx _ 



a2 q > q 



X 



Both versions of this bi-fractal model require three parameters. The multi- 
fractal spectrum is degenerated to two single points, thus its width can be 
defined as Aa = ai — a2- 
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9 Future Directions 



The most straightforward future direction is to analyse more types of time 
series from other complex systems than those listed in Chapter 2 to check 
for the presence of fractal scaling and in particular long-term correlations. 
Such applications may include (i) data that are not recorded as function of 
time but as function of another parameter and (ii) higher dimensional data. 
In particular, the inter-relationship between fractal time series and spatially 
fractal structures can be studied. Studies of fields with fractal scaling in 
time and space have already been performed in Geophysics. In some cases 
studying new types of data will require dealing with more difficult types of 
non-stationarities and transient behaviour, making further development of 
the methods necessary. In many studies, detrending methods have not been 
applied yet. However, discovering fractal scaling in more and more systems 
cannot be an aim on its own. 

Up to now, the reasons for observed fractal or multifractal scaling are 
not clear in most applications. It is thus highly desirable to study causes 
for fractal and multifractal correlations in time series, which is a difficult 
task, of course. One approach might be based on modelling and comparing 
the fractal aspects of real and modelled time series by applying the methods 
described in this article. The fractal or multifractal characterisation can 
thus be helpful in improving the models. For many applications, practically 
usable models which display fractal or transient fractal scaling still have to 
be developed. One example for a model explaining fractal scaling might 
be a precipitation, storage, and runoff model, in which the fractal scaling 
of runoff time series could be explained by fractional integration of rainfall 
in soil, groundwater reservoirs, or river networks characterised by a fractal 
structure. Also studies regarding the inter-relation between fractal scaling 
and complex networks, representing the structure of a complex system, are 
desirable. This way one could gain an interpretation of the causes for fractal 
behaviour. 

Another direction of future research is regarding the linear and especially 
non-linear inter-relations between several time series. There is great need 
for improved methods characterising cross-correlations and similar statistical 
inter-relations between several non- stationary time series. Most methods 
available so far are reserved to stationary data, which is, however, hardly 
found in natural recordings. See |101j for a very recent approach analysing 
fractal cross-correlations in non- stationary data. An even more ambitious 
aim is the (time-dependent) characterisation of a larger network of signals. 
In such a network, the signals themselves would represent the nodes, while 
the (possibly directed) inter-relations between each pair represent the links 
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(or bonds) between the nodes. The properties of both, nodes and hnks can 
vary with time or change abruptly, when the represented complex system 
goes through a phase transition. 

Finally, more work will have to be invested in studying the practical 
consequences of fractal scaling in time series. Studies should particularly 
focus on predictions of future values and behaviour of time series and whole 
complex systems. This is very relevant, not only in hydrology and climate 
research, where a clear distinguishing of trends and natural fluctuations is 
crucial, but also for predicting dangerous medical events on-line in patients 
based on the continuous recording of time series. 
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Figure 1: (a) Examples of self-affine series Xi characterized by different Hurst exponents H — 0.9, 0.7, 0.5, 0.3 (from top to bottom). The data has 
been generated by Fourier fiUering (see Section 8. 1) using the same seed for the random number generator, (b) Differentiated series Axi of the data 
from (a); the Ax, are chai'acterized by positive long-term con'elations (persistence) with y — 0.2 and 0.6 (top curve and second curve), uncoiTelated 
behaviour (third curve), and anti-correlations (bottom curve), respectively. 
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Figure 2: Compaiison of the autoconelation functions C(s) (decreasing functions) and fluctuation functions Fjis) (increasing functions) for short- 
term correlated data (top) and long-term correlated data (y = 0.4, bottom). The asymptotic slope H~a = 0.5 of Fiis) clearly indicates the absence of 
long-term correlations, while H~a—\—y/2> 0.5 indicates the presence long-term correlations. The difference is much hai'der to observe in C(s), 
where there are more statistical fluctuations and even negative values (e.g., above s - 150 in (a) and between s — 300 and 400 in (b), not shown). 
The data have been generated by an AR process (Eq. (4)) and Fourier filtering (Section 8. 1 ) for (a) and (b), respectively. The dashed lines indicate 
the theoretical curves. 
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Figure 3: Spectral analysis of a fractal time series characterized by long-term conelations with y = 0.4 (P = 0.6). The expected scaling behaviour 
(dashed line indicating the slope -[3) is observed only after binning of the spectrum (circles). The data has been generated by Fourier filtering (see 
Section 8.1). 
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• Uncorreiated Gaussian data (20 systems) 

■ LTC Gaussian data (h^ = 0.6) (20 systems) 

♦ Sum of 3 AR processes with correlation times 3, 21, and 500 - 
A Extreme configuration of the sum of AR processes (out of 20) 
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Figure 5: Application of discrete wavelet transform (WT) analysis on uncorreiated data (black circles), long-term correlated data (y- 0.8, a - 0.6, 
red squares), and short-term correlated data (summation of three AR processes, green diamonds). Averages of Fiis) averaged over 20 series with A' 
- 2^^ points and divided by s^'~ are shown, so that a horizontal line corresponds to uncorreiated behaviour. The blue open triangles show the result 
for one selected extreme configuration, where it is hard to decide about the existence of long-term correlations (figure prepared by Mikhail 
Bogachev). 
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Figure 7: Example of the Wavelet Transform Modulus Maxima (WTMM) method, showing the original data (top), its eontinuous wavelet 
transform (grey scale coded amplitude of wavelet coefficients, middle), and the extracted maxima lines (bottom) (figure taken from [67]). 
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Figure 8: Multifractal Detrended Fluctuation Analysis (MF-DFA) of data from the binomial multifractal model (see Section 8.3) with a = 0.75. 
F,i{s) is plotted versus s for the q values given in the legend; the slopes of the cui'ves correspond to the values of h{q). The dashed lines have the 
slopes of the theoretical slopes h(±oa) from Eq. (42). 100 configurations have been averaged. 



56 



1± 



10 



20 



30 



Figure 9: Illustration for the definition of return intervals between extreme events above two quantiles (thresholds) q\ and q2 (figure by Jan 
Eichner). 
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Figure 10: Normalized rescaled distribution density functions R,,P,/(r) of r values with Re/ = 100 as a function of r/R,/ for long-teim correlated data 
with Y- 0.4 (open symbols) and y- 0.2 (filled symbols; we multiplied the data for the filled symbols by a factor 100 to avoid overlapping cui'ves). 
In (a) the original data was Gaussian distributed, in (b) exponentially distributed, in (c) power-law distributed with power -5.5, and in (d) log- 
normally distributed. All four figures follow quite well stretched exponential curves (solid lines) over several decades. For small r/R^ values a 
power-law regime seems to dominate, while on large scales deviations from the stretched exponential behaviour are due to finite-size effects (figure 
by Jan Eichner). 
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Figure 11: Illustration for the definition of maxima niR within periods of R- 365 values (figure by Jan Eichner). 
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